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Periodically driven quantum systems can be nsed to realize quantnm pumps, ratchets, artificial 
gauge fields and novel topological states of matter. Starting from the Keldysh approach, we de¬ 
velop a formalism, the Floquet-Boltzmann equation, to describe the dynamics and the scattering of 
quasiparticles in such systems. The theory builds on a separation of time-scales. Rapid, periodic 
oscillations occurring on a time scale To = 27r/D, are treated using the Floquet formalism and 
quasiparticles are defined as eigenstates of a non-interacting Floquet Hamiltonian. The dynamics 
on much longer time scales, however, is modelled by a Boltzmann equation which describes the 
semiclassical dynamics of the Floquet-quasiparticles and their scattering processes. As the energy 
is conserved only modulo Tiff, the interacting system heats up in the long-time limit. As a hrst ap¬ 
plication of this approach, we compute the heating rate for a cold-atom system, where a periodical 
shaking of the lattice was used to realize the Haldane model pQ . 


Periodically modulated quantum systems can effec¬ 
tively be described by a static Hamiltonian. This the¬ 
oretical concept has recently evolved into a major exper¬ 
imental tool used by many groups to generate new states 
of matter. 

Early experiments ElE] used, for example, that one 
can effectively change the strength as well as sign of the 
hopping of atoms in an optical lattice, allowing to realize 
new types of band structures. Periodic driving has also 
be used to realize directed transport in quantum ratch¬ 
ets [4] . More recently, the realization of emergent Gauge 
fields and topological band structures has been at the 
focus of many studies. Examples of such experiments 
include the generation of Gauge fields and superffuids 
with finite momentum[5l[6], the generation of topological 
quantum walks |7] and of effective electric fields in a dis¬ 
crete quantum simulator [S] , the realization of (Floquet-) 
topological insulators with photons [9] , the generation of 
spin-orbit coupling cni, the direct measurement of Ghern 
numbers and Berry phases in the Hofstadter Hamiltonian 
[mill] and the realization of quantum pumps [ig. A 
recent experiment of the Esslinger group beautifully re¬ 
alized the Haldane model [T], i.e., a model which demon¬ 
strates that a quantum-Hall state can exist without ho¬ 
mogeneous external magnetic fields. In solid-state sys¬ 
tems circularly polarized light has been used |14j to ma¬ 
nipulate the surface states of topological insulators. 

Also from the theory side, many proposals have 
pointed out that periodically driven states can be used to 
realized a wide range of states of matter. Examples are 
photoinduced quantum Hall states [miis] and various 
topological Floquet states [IM] including dissipative 
systems [21], quantum ratchets for Mott insulators [22] . 
Majorana Fermions in driven quantum wires [23], the 
generation of non-abelian gauge fields [23] , Floquet frac¬ 
tional Chern insulators [2S] , Floquet-Anderson insulators 
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and quantized charge pumps 

In a periodically driven system, the Hamiltonian has 
only a discrete time-translational symmetry, H{t + TQ) = 
H {t). As a consequence, the total energy is not conserved 
but quantized changes of energy are possible, AE = n hfl 
with H = 2tt/To and n € Wi. For non-interacting sys¬ 
tems the absence of energy conservation has mostly no 
effect. The situation is, however, different when inter¬ 
actions in a many-particle system are considered. For a 
generic closed system, one can expect that in the long¬ 
time limit, t —>■ oo, the system approaches the state 
with the highest entropy consistent with the conserva¬ 
tion laws. In the absence of some cooling mechanism, 
e.g., by an external bath m or by emitting radiation, one 
can therefore expect that generic interacting Floquet sys¬ 
tems heat up to infinite temperatures [28] (an exception 
are many-body localized systems [29]). This important 
(and well-known) aspect has received relatively little at¬ 
tention in previous studies. Eckstein and Werner used 
time-dependent dynamical mean-field theory to study 
heating effects. In m, the stability of BEG conden¬ 
sates in periodically driven systems discussed based on 
phase-space arguments for energy non-conserving scat¬ 
tering processes. 

The goal of this paper is to derive and apply a Floquet- 
Boltzmann equation, i.e., a kinetic equation which can 
be used to describe the dynamics of weakly interacting 
Floquet systems. Such a kinetic equation is perhaps 
the simplest theoretical description which captures mi¬ 
croscopically how interactions can equilibrate an inter¬ 
acting quantum system. Like other quasi-classical ki¬ 
netic equations, our approach builds on a separation of 
time scales and describes situations, where the change 
of occupation functions is much slower than Tq and the 
time-scales set by the bare parameters of the Hamilto¬ 
nian. Motivated by the recent realization of the Haldane 
model by the Esslinger group [1], we investigate the ef¬ 
fects of local interactions in a fermionic system. Starting 
from the Keldysh-formalism, we obtain a quantum ki- 



2 


netic equation from the Dyson equation, which is then 
reduced to the Floquet-Boltzmann equation. As an ex¬ 
ample, we quantitatively investigate heating rates in a 
limit when energy-conserving processes dominate (real¬ 
ized in Ref. ID)- 

A similar semiclassical kinetic equation was also de¬ 
rived in a recent preprint of Seetharam et al. [27) . In 
contrast to our work, they considered electron-phonon 
coupling instead of fermion-fermion interactions. Tech¬ 
nically, Ref. m used an equation-of-motion approach. 
While this approach is, perhaps, less well suited to in¬ 
vestigate limitations of semiclassics, we expect that in 
the semiclassical limit it gives results equivalent to our 
derivation. 


I. KELDYSH APPROACH AND QUANTUM 
KINETIC EQUATION 

Our aim is to develop a Floquet-Boltzmann approach 
for interacting many-particle systems. The derivation 
builds on two main elements: a separation of time scales 
and the limit of weak interactions. We consider a time- 
dependent many-particle Hamiltonian H{t) which has 
the property that on short, microscopic time scales it 
is approximately periodic, 

Hit)^H{t + To) ( 1 ) 

with period Tq. Furthermore, we also allow for a slow 
time-dependence on time-scales large compared to the 
To and all microscopic time scales (inverse kinetic and 
potential energies and inverse scattering rates). The lat¬ 
ter can, for example, be used to describe the influence of 
external forces or the slow (quasi-adiabatic) change of the 
Hamiltonian. The Hamiltonian can therefore be written 
in the form 

H{t) = (2) 

n 

with D = ^ and with Fourier series coefficients (t) = 
H-n{ty which are time-dependent on time scales much 
larger than Tq. 

As the derivation of quantum kinetic equations is ulti¬ 
mately based on perturbation theory, we further require 
that this perturbation theory can indeed be applied. In 
our examples, this will be the case when interactions are 
weak. More generally, one can use the approach also in 
cases where interactions are strong, but scattering rates 
are nevertheless weak either due to phase space restric¬ 
tion (e.g., for the fermionic quasiparticles of a Fermi liq¬ 
uid close to the Fermi surface) or just because the density 
of excitations is low (e.g., a weakly excited bosonic Mott 
insulator [22)1. 

To simplify the presentation, we will not discuss the 
most general setup but restrict ourselves to a simpler 
case. We consider weakly interacting Fermions in a lat¬ 
tice model that is described by the following Hamiltonian 

= H0(t)(3) 



FIG. 1. (Color online) Two of four interaction vertices gen¬ 
erated by interaction of type Eq. Q. The remaining two 
diagrams are obtained by complex conjugation of (a) and (b), 
which simply corresponds to inverting the direction of the 
arrows. 


with a static and local interaction 

= uj2 
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We assume in the following that the interaction strength 
U is sufficiently small, i.e., 17 ^ J, in order to allow a per¬ 
turbative approach to solve the problem. Note that our 
approach can be generalized in a straightforward way to 
more complicated and time-dependent interactions. We 
will present a theory that will unite the Floquet theory 
describing the oscillatory character of our model with 
the Keldysh approach of quantum field theories captur¬ 
ing the non-equilibrium behaviour of the system due to 
interactions and adiabatic drifts. 

In order to derive a quantum kinetic equation, we use 
the standard Keldysh approach [3T| [32]. We will not 
give a review of this Keldysh approach here, but just 
give a few main definitions. More details can, e.g., be 
found in the book by Kamenev [31] on which the following 
discussion is based. 

To describe the time-evolution of the density matrix, 
p{t) = U{t)poU{ty one needs to keep track of two 
time-evolution operators U{t) and U{t)^. Within the 
functional-integral version of the Keldysh approach one 
therefore introduces Grassmann fields '0 on the forward 
branch, , and on the backward branch, 'tp~ . One then 
performs a rotation in Keldysh space using the following 
relations 



i k] 
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( 5 ) 

After such a Keldysh rotation the action corresponding 
to the system excluding the interactions, i.e., H — i7i„t, 
can be written in the form 


^0 = 


/ 





( 6 ) 


where the integral should be understood as / = 

J dtj dt' and the fields as = ipi{t) and ijj = 

The non-interacting Green’s function, G, is a 2x2 matrix 
with components 


Go' = 


(Go')'' (Go"')'" 


(G°o')^ 


( 7 ) 
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Generally, the Green’s functions of the full system are 
given by 

G^{x,x') = -i0{t - t'){{tfj{x),tfj'^{x')}) 
G^{x,x') = ie{t' - t)({V’(a;),(a;')}) 

G^{x, x') = -i{[^{x), ik^x')]) (8) 

with X = {r, t) or within the functional integral by 

G°‘^{x,x') = -i [ (9) 


Note that we choose ?i = 1 throughout the paper. Here, 
G^^ = G^ and G^^ = G^ are the retarded and advanced 
Green’s functions, where G^ = (G^)^, and G^^ = G^ is 
the Keldysh Green’s function (note that G^^ =0). The 
latter can generally be parametrized in the following form 

G^ = G^oF-FoG^ (10) 

where F = F{x, x') is a Hermitian matrix and called 
distribution matrix. The “o” represents matrix multipli¬ 
cation in all indices (space, time, spin). While G^,G^ 
carry information about the spectrum of the system, F 
holds the information about the occupation as we will 
see below. 

The Dyson equation 

(Gq-i - S[G]) o G = i (11) 


plays a central role in the derivation of the quantum ki¬ 
netic equation. All objects are again matrices in Keldysh 
space, and the “o” now also includes a matrix multipli¬ 
cation with respect to Keldysh indices. Here, one uses 
that the self-energy, S = S[G], can be viewed as a func¬ 
tional of the full Green’s function, G. The Dyson equa¬ 
tion therefore is an integro-differental equation to deter¬ 
mine G in a time dependent system. The Keldysh com¬ 
ponent of the equation above can be identified with the 
quantum kinetic equation (see below). Upon further ap¬ 
proximations, this equation can be simplified to obtain 
the semiclassical Boltzmann equation. 

Using the fact that (in the fermionic case) the self¬ 
energy has the same structure as G~^ and G, i.e.. 


E = 




( 12 ) 


one can write down the Keldysh component of the Dyson 
equation, i.e. the quantum kinetic equation, 

Fo{G^)-^-{G^)-^oF = E^-(S^oF-FoE'^) (13) 

So far the expression is exact and no approximation has 
been made. For an interaction of the type as in Eq.Q, 
the interaction part of the action takes in Keldysh space 
the form 

'S'int =-^ j (14) 

2,(T _ _ 



FIG. 2. (Color online) Second order diagrams due to on-site 
interactions, (a)-(e) contribute to and (f)-(l) contribute 
to E^. E^ is obtained by realising that E'^ = E^F Note 
that the assignment of the spin indices as in (a) is the same 
for all diagrams. 


where a =l (f) for tr =t (i). Diagrammatically, 
S'int yields four interaction vertices: two independent 
ones plus their complex conjugates (see FigQ. To de¬ 
rive a quantum kinetic equation, we will consider (self- 
consistent) self-energy diagrams up to second order. 

To linear order in U, one obtains the familiar Hartree- 
Fock contributions: the energy of an f electron is changed 
by U{nii). In the following, we will absorb all these com¬ 
pletely into a redefinition of the non-interacting part 
Note, however, that due to the time dependence of expec¬ 
tation values (nil), the non-interacting part will obtain 
an extra time-dependence which has to be calculated self- 
consistently. 

Due to the fact that G^^ = 0, twelve indepen¬ 
dent diagrams are obtained in total for the sec¬ 
ond order expansion (see Fig§ that contribute 
to the respective part of the self-energy. Ex¬ 
ploiting the general identities for retarded and 
advanced Green’s functions G^{x', x)G^{x, x') + 
G^{x’ ,x)G^{x,x’) = —G^^{x’,x)G^"^{x,x’) and 

{G^{x,x')Y + {G^{x,x')Y = {G^'^{x,x’)Y, where 

GRA ^ qR _ qA^ 

one can write down the second order 
contributions to the individual parts of the self-energy 
as 


E 


K 

a 


- G^{x',x) G^{x,x') G^{x,x') (15) 

+ G’^{x',x) G^^{x,x') G^^{x,x') 

- G^^{x',x) G^{x,x') G^^{x,x') 

- G§^{x',x) G^^{x,x') G’^{x,x') 
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S? = - 


U' 


2 r 


(16) 


^ ^Gi{x',x)G^^{x,x')G^^{x,x') 

+ Gi{x',x) G^{x,x') G^{x,x') 

+ G^{x',x) G^{x,x') G^{x,x') 

+ G^{x\x) G^(x,x') G^{x,x') 


where again tr =4, (f) for a (1) is the conjugate spin 
to cr and T,^ = Yi^^{x, x'). The expression for can be 
straightforwardly obtained by using that 

The equations for the self-energies as func¬ 

tions of the interacting Green’s function together with 
the Dyson equation (11) and a suitable initial condition 


define a quantum kinetic equation (QKE), which can be 
used to study the dynamics of a time-dependent inter¬ 
acting system. 

A direct (numerical) solution of the QKE is very 
challenging, especially for a time-dependent Hamilto¬ 
nian where the Green’s function depends on two time- 
variables separately. One can, however, make progress in 
situations where there is a clear separation of time scales. 

The next step is to switch to a semi-classical represen¬ 
tation via the Wigner transformation which we have to 
combine with the Eloquet formalism to take into account 
rapid periodic oscillations. 


II. FLOQUET EIGENSTATES AND ELOQUET 
GREEN’S FUNGTIONS 

The analysis of the dynamics of the Eloquet system 
starts from the non-interacting, but time-dependent part 
of the Hamiltonian, H°{t) = Y. , where 

iJ°(t) varies only on time scales Tgiow 3> Tq. More pre¬ 
cisely, we will include in also all Hartree-Eock cor¬ 
rections arising from the interacting part of the Hamilto¬ 
nian, i.e., terms like 17 

where {rLi„{t)) will generically have components oscil¬ 
lating with frequencies nVL plus an extra slow time- 
dependence arising, e.g., from heating. Using the as¬ 
sumed separation of time scales allows to define 

<(t) = ^«(to)),,4c^.e-^"* (17) 

as an approximation to H^{t) for t close to to (on the time 
scale set by TsIow) with + Tq) = To solve 

such exact time-periodic Hamiltonians one uses the Elo¬ 
quet theorem (the analog of the Bloch theorem for peri¬ 
odic time- instead of space-dependencies) stating that so¬ 
lutions can be written in the form \4’v{t)) = e“*'’‘'*|(/)j,(t)). 
Here, \4’iy{i)) = \4’v{t + 7b)) is a time-periodic Eloquet 
state and is called the quasi energy. The eigenstates 
to) are obtained by diagonalizing the Eloquet Hamil¬ 
tonian 

^to,nm ~ 'IT'^Snm 1 + ( 18 ) 

with to). Eor practical calculations 

, the Eloquet indices n,m run from —Nf to Nf, which 


is chosen such that flNf is much larger than any other 
energy scale in the problem. Due to the ‘translation’ 
invariance (obtained for Nf —>■ oo), = D -|- 

rmi on.® ®nn obtain from each eigenstate with energy e 
an eigenstate with energy by a simple translation of 

the Eloquet indices, n ^ n + k. Therefore, it is sufficient 
to consider only eigenstates with eigenenergies 


n n 

~ 2 - ^*0'' ^ 2 


(19) 


Here, ly encodes the usual quantum number (spin, band- 
index, momentum). 

The Heisenberg operator 




4 ( 20 ) 


creates a fermion in such a Eloquet eigenstate with 
!^(*) = (here i includes all quantum num¬ 

bers, e.g., lattice site and spin). Note that 

ifLMftoA*')) ( 21 ) 

the occupation of the Eloquet states, is time indepen¬ 
dent for the non-interacting, periodic Eloquet Hamilto¬ 
nian ( [T7| ). 

Here, it is important to stress that a single function 
y describes the occupation of the Eloquet states, and 
it is not necessary to introduce separate occupations for 
each Eloquet index (further Eloquet eigenstates obtained 
by translations in Eloquet space yield exactly the same 
Eloquet-creation operator, j/(t))- Our main goal will 
be to find a semiclassical description of the time evolution 
of the Eloquet-occupations nt^u{t). 

In a state with {fl„^^{t)f^^y{t)) = one can 

calculate the Green’s functions (|^ of the noninteracting 
system I©- They depend on two time variables but the 
dependence on {t-\-t')/2 is purely periodic and therefore 
it is useful to represent the Green’s function in Eloquet 
space by introducing the Eloquet-Eourier transformation 
Anm{io) = 4/ to ob¬ 

tain 


(^R . ^ <(’"o!)^(^) 

nm ^ W - etg 1/ - ZD-I-lO 

y 

y w - fftg ^ - ZD - lO 
^^.b(‘*^) = “ 27 rzy b(w- ZD) x 


nm 




v,l 


X <by(*)C:t'(j)(l-2rqg..) 


j.m+1 


( 22 ) 


The distribution function defined in ( |10[ ) is therefore 
given in Eloquet representation by 


FtoM = E C>^(*) C:t'(j) (1 - 2n,g..) (23) 

nm 

u,l 
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Note that the Green’s functions defined above are solu¬ 
tions of the Dyson equation 

'y ~ ^to,nm')^to,m'm{^,t) = Snm^ (24) 

m' 

III. FLOQUET-WIGNER FORMALISM AND 

FLOQUET-MOYAL EXPANSION 

Our central goal is to separate the slow dynamics, 
treated within a semiclassical approach, from rapid, peri¬ 
odic oscillations which have to be treated fully quantum- 
mechanically. The latter aspect is treated within the Flo- 
quet formalism (see Sec. 0 which uses that in a strictly 
periodic system with period Tq, only Fourier modes of 
the form nD with D = occur. To derive the slow, 
semiclassical dynamics, the starting point is the use of 
a Wigner representation of the Green’s function 
usually obtained by using a Fourier transformation of the 
relative time coordinate, t — t', for fixed t = 

These two approaches can be combined in situations 
where there is a clear separation of time scales, as dis¬ 
cussed above. We require that the oscillation period Tq 
is much smaller than all time scales, Tgiow, on which the 
occupation function changes or on which the oscillating 
Hamiltonian is modified. This allows us to introduce the 
following ‘Floquet-Wigner representation’ for functions 
A{t,f) (Green’s functions or self-energies) which depend 
on two time coordinates 

— (2fi) 

Here, 

<5#- t) = —(26) 

VSttt 

is a version of the i5-function which is broadened on the 
time scale r chosen much larger than the period Tq, but 
much smaller than the time scale of slow modifications, 

'^slow 

To < T < Tsiow (27) 

The use of the filter function guarantees 

that is not rapidly oscillating on the time 

scale Tq. More precisely, all oscillating components at 
frequency D are exponentially suppressed by the factor 
g-(f2'r)’^/4 _ {t/T o)^ ^]pg convolution with the 

filter function. Later, i will take over the role of to intro¬ 
duced in the previous section. 

The inverse transformation of the Floquet-Wigner rep¬ 
resentation is given by 

A{t,t')^Yl / ^ A„„(cc,^) 

n^m 

(28) 


Due to the finite width r of the filter function (^^(t), this 
back-transformation is not exact but it is valid (with ex¬ 
ponential precision) in situations when Eq. ( |27[) h olds. 
To see just this, it is instructive to plug Eq. (|28p into 
Eq. (251: The condition Tq guarantees that the Flo- 


quet indices do not mix, while r ^ TsIow allows to use 
use 5r(^) as a true J-function for all components which 
vary slowly in time. 

By definition, the u argument of the Floquet-Wigner 
representation is restricted to the interval — § < w < ^ 
(the analog of the reduced Brillouin zone for periodic 
systems in real space). 


Within the quantum kinetic equation (131, one has to 
compute the product of matrices, A = BoC, which takes 
the form 

= [B o C){t2,ti) = J dt'B{t2,t')C{t',ti) (29) 


Into this equation we plug the inverse Floquet-Wigner 
transformation, Eq. (28), for B and C, and Taylor- 
expand the time and frequency arguments of Bnmi^A) 
and t). It turns out, that the resulting expression 

can be written in the form 


(30) 

Here denotes the time- or frequency derivative of 

the functions A or _B, respectively. By a Taylor-expansion 
of the exponential, one obtains the well-known Moyal 
expansion m- The only difference in comparison to the 
standard Moyal expansion is that it is supplemented by 
a simple matrix multiplication of the Floquet indices. A 
short derivation of the Floquet-Moyal expansion can be 
found in appendix [A} 

For the computation of the self-energy, we also need 
the Floquet-Wigner transformation of a different type of 
product given by A{t,t') = B{t,t')C{t,t'). In this case 
one obtains directly 

/ duj^ 

Bn'm'{u}',t) (31) 

X Cn—n',ra—m' (^ ^ A) 


IV. FLOQUET-BOLTZMANN EQUATION 

Boltzmann equations are a powerful tool to describe 
how scattering processes affect the semiclassical dynam¬ 
ics. They do not aim at describing quantum-coherent 
processes at short times, but instead focus on the physics 
at time scales set by slow changes of external parame¬ 
ters and by the scattering time of particles. It there¬ 
fore builds on a clear separation of the time scales for 
quantum-coherent processes (captured by us within the 
Floquet approach for periodically driven system) and 
for the semiclassical dynamics which changes occupation 
functions. 
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The derivation of the Floquet-Boltzmann equation can 
be divided into four steps, (i) Starting point is the quan¬ 
tum kinetic equation (Hi together with the calculation of 
self-energy diagrams, which are functionals of the Green’s 
function, see Eq. The next goal is to use the 

separation of time scales. Therefore, (ii) one uses the 
Floquet-Wigner representation, introduced in Sec. m 
for all Green’s functions and self energies. Gonvolutions, 
’o’, can be written in terms of a Floquet-Moyal product, 
(30). Using the separation of time scales which implies 
that terms proportional to dt give small contributions, 
we can (hi) Taylor-expand ) to leading 

order, i.e., to linear order on the left-hand and to ze¬ 
roth order on the right-hand side of the quantum kinetic 
equation (13). For problems which are not spatially ho¬ 
mogeneous a similar Moyal expansion is also used for the 
spatial coordinates. Finally, (iv) the resulting equation is 
projected onto on-shell processes, e.g., by an integration 
over frequencies. 

The final result of these steps is an equation for the oc¬ 
cupation functions nk,^ (r, t) of the Floquet eigenstates at 
time t. Here k is the momentum and ^ includes band- and 
spin indices. The Floquet states at time t are the eigen¬ 
states, of the Floquet Hamiltonian (18) with 

1 / = (k, and we have to set to = t- In complete analogy 
to the treatment of the variable to = t-, which we used 
in Sec. [H] to deal with the slow time dependence, we also 
allow that the Hamiltonian depends smoothly on the spa¬ 
tial parameter Tq = r. For a lattice model with sites 
per unit cell, the eigenfunctions are calculated in momen¬ 
tum space by diagonalizing a n„(2iV/ -I- 1) x n„(2iV/ -|-1) 
dimensional matrix. In the following we will denote 
the corresponding eigenfunctions in momentum space by 
^tr k {(0 where n = —Nf,—Nf -|- 1 ,..., Nf is the Floquet 
index and i = l,...,nu describes the structure of the 
Bloch-Floquet wave function within the unit cell. Some¬ 
times, we will omit the t and r index to simplify notations 
and just write ^{i). 

For the following discussion, we will not discuss the 
(main) part of the derivation which is identical for Flo¬ 
quet systems and conventional cases, as these are well 
described in the literature [311 133] and textbooks m- 
Instead, we will only describe those aspects which are 
different in the Floquet case. 


A. Semiclassical dynamics and left-hand side of the 
Floquet-Boltzmann equation 


The Floquet-Moyal expansion, Eq. (30), differs from 


the standard Moyal expansion only by the presence of 
the extra Floquet indices. This gives rise to a simple ma¬ 
trix multiplication. To be able to describe also situations 
where the occupation functions are not spatially trans¬ 
lational invariant, but vary smoothly (on length scales 
large compared to the lattice spacing), one uses a Wigner 
representation and Moyal expansion for space and mo¬ 
mentum degrees of freedom, similar to the one described 


above for time and frequency variables [31H33j . A major 
difference between the momentum and the frequency de¬ 
pendence is, however, that the semiclassical occupation 
functions depend on the quantum numbers momentum 
and band index, but not on frequency and Floquet in¬ 
dices. 

A derivation of the conventional collisionless Boltz¬ 
mann equation based on the quantum kinetic approach, 
which (in contrast to previous derivations) includes all 
Berry-phase correction to leading-order, has recently 
been given by Wickels and Belzig [33]. One can check 
(see Appendix]^ that the only difference arising in the 
Floquet case is that all Berry curvatures have to be 
computed from the Floquet-eigenfunctions introduced in 
Sec. [Hi 

(32) 

i,n 

where the scalar product involves also a summation 
over the Floquet index n. Here 9^ and 9^, /r, iz = 
(t,T’l)■'’ 2 ,?' 3 ,Pi,P 2 ,P 3 ) stands for derivatives in time, 
space and momentum variables |34j . These Berry cur¬ 
vatures modify the semiclassical equations of motion in 
phase space and therefore also the left-hand side of the 
Boltzmann equation which takes the form 

{dt T^k.cVfc -f Uk,5Vr)nk,{(r,t) = 9tnk,{(t)|^„ii (33) 
with [331133] 

Uk,{ = (1 -|- Jl^p) • Vk£{ — + ^pp ■ Vr£^ 

■^k,{ = ~(1 + ^rp) ■ + ^rr ' ^k£4 (34) 

While rJrp, and Opp are matrices with r,p = (1, 2,3), 
Clj-t and ttpt are vectors {t = 1). Note also that ftm ^pp 
can be related to effective magnetic fields, and ^pt 
are referred to effective electric fields. 

Since many modern applications of Floquet Hamilto¬ 
nians [T] E] EIHIl [13] have as a goal to realize systems 
with non-trivial Berry phases, it is important to keep 
track of these effects on the left-hand side of the Boltz¬ 
mann equation. Gonsider, for example, an interacting 
Floquet system which heats up as function of time (see 
Sec. |V]). The (slow) change of occupation functions can 
trigger a change of the momentum-space Berry curvature 
r2piP2, a momentum-space ‘magnetic’ field. This implies 
that also corresponding momentum-space ‘electric’ fields 
fltpi and H(p 2 are generated. They can, e.g., induce a 
macroscopic rotation of the cold-atom system. 


B. Scattering and the right-hand side of the 
Floquet-Boltzmann equation 


To calculate the right-hand side of the Floquet- 
Boltzmann equation, we start from the self-energies 
Eq. (15) and (16). First, we need an expression for the 
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Green’s function. Due to the assumed separation of time 
scales, it is sufficient to evaluate the Green’s function us¬ 
ing a zeroth-order Floquet-Moyal expansion of the Dyson 
equation (111. Furthermore, within our perturbative ap¬ 


proach, we do not have to include any self-energy cor¬ 
rections (as the self-energy is already (x U^). Using the 
Floquet-Wigner representation of both the Hamiltonian 
and the Green’s function, the Dyson equation takes with 
these approximations exactly the form of Eq. (24|. This 


implies that we are allowed to use directly the Green’s 
functions of Eq. (22) with n^, = Uk ^(t) and etg^u = eqk.^- 
To evaluate the Eloquet-Wigner representation of the 
self-energies Eq. (15),(16l and the right-hand side of 
quantum-kinetic equation (13), we use the convolution 


formula Eq. (31) twice. For the first line of the formula 


in Eq. (15), for example, we obtain after a few steps of 


simplihcation a contribution of the form 


^K,l 

'ij,C7 
nm 

iirU'^ 


(35) 


2 ^ ^ j {2n/aY {2n/aY 

X S(^UJ -|- Cq—^k--p,?7,CT ^q,A,fT “1“ ^ 'u)r^) 

^ (1 27iq_p_^ ct )(1 ^TZk— p,? 7 ,ct ) (1 2 ?lq , 7 ) 

where a =| (t) for a =t (i), the superscript ’1’ refers to 
the first line of Eq. ( [I^ , and 


kpq,t 

im' -\-s 


m' ,m" 
Trn^' -\-s 


-ij) 

k— p,?7,<tv / ^k—p,?7,<7 / 




(36) 


ample, the form 
f duj 


27r 


Tr{Hk.5(a;)E^'i(a;,k)} 


(37) 


r;,/ 2 ,A l,s,u ‘ 


dp 


dq 


(27r/a)'^ (27r/a)‘’* 




p,/i,(7 ^k—p,77,fT ^q,A,(T “b ^s/uU) 

X (1 2Uq_p^^ cr') (^ p.r^.crO ^^q.AjCr) 

with Asiu = s—l—u, and the transformed matrix element 

$ 


iAA^) 

= E 


kpq 

ij,nm 

kpq 



(38) 


It is convenient to introduce for each occupation function 
a separate momentum variable and a d-function which 
guarantees momentum conservation (modulo reciprocal 
lattice vectors Ga). Performing the entire procedure for 
all terms of Eq. (15) and likewise for all contributions 


associated with the second term on the right-hand side of 
the quantum kinetic equation (13), one eventually hnds 


an expression for the collision integral 


Aoiil^k.j.o-] — 

77 ,/i,A oc,n 

X WP 


dqi dqa dqg 
{2tt/ aY {2111aY {2iT/aY 


jprjA.a (27r/a) d(k -b qi - qa - qa - aG) 

kqiq2q3 


X ^(ek,^,^ + 




— e, 


q2,r;,(T £q3,A,o- 




q2,7?,^7’ 


^q3,A,(T (1 ^k,^,cr) (1 


— nD) 

qi.M.ff) 


^k,^,fT ^qi,/i,CT (1 


^q2,i;,o-) (1 ’Aqj^^^o-)] 

(39) 


where /, s, u, n, m, m', m" are Eloquet indices, rj, A are 
band indices and i,j denote sites within the unit cell. 
Here, {2iTjaY is the volume of the Brillouin zone. We 
have omitted extra r and t labels which each function 
obtains to reflect the smooth time and spatial dependen¬ 
cies of the system. 

The last remaining step is to evaluate the resulting 
formula on-shell: we multiply the right-hand side of the 
quantum kinetic equation (13) by the Floquet-spectral 
function 24k,{,nm(w) « 2-k5(Yj~ ^he 

state with quantum numbers k and integrate over fre¬ 
quencies and trace over Floquet- and space indices. Note 
that considering only diagonal, on-shell contributions im¬ 
plies that Boltzmann-type equations cannot describe co¬ 
herent quantum-oscillations. The resulting equations are 
therefore only valid on time scales longer than the de¬ 
cay time of such oscillations. This is consistent with our 
assumptions on the separation of time scales underly¬ 
ing our analysis. Furthermore, a quasiparticle has to be 
well defined, implying that the broadening of the spectral 
function by scattering is small compared to the energy of 
the quasiparticles (and therefore also small compared to 
D). 


where we have introduced the integers a,n G S to ac¬ 
count for Umklapp scattering in momentum- and fre¬ 
quency space, respectively, and ^ is the scatter- 

kqiq2q3 

ing rate for a process involving an energy transfer to the 
system of nil, n GTL. We obtain 




kqiq2qs 


= 2ttU^ 




I«qiq2q3 


(40) 


—( 711+712 — 713 — 724 ) 


(41) 


with the amplitude 

rii., = E 

^qiq 2 q 3 1,711,712,713,714 

X C?.x*) KUA^) AAA^) 

Note that the Floquet- and momentum indices enter the 
matrix elements, and therefore the collision integral, in 
a completely different way: occupation functions depend 
on momentum and band indices, but do not depend on 
the Floquet indices. Correspondingly, we sum over Flo¬ 
quet indices in Eq. (411, but not over momentum or band 


After this last transformation Eq. (35) takes, for ex¬ 


indices. We will discuss this important difference again 
in the concluding section. 



































The collision integral in Eq. (391 forms the right-hand 


side of the Floquet-Boltzmann equation 

F T )^k,^,cr (1*5 — 'f^coll[^k,^,a '(^7 ^)] 

(42) 

where, in general, also the effective forces and velocities 
depend smoothly on time and space, = ^k,{(r,i) 
and Vk.j = i’k,^(r, t). This dependence can either arise 
from an explicit r and t dependence of the Hamiltonian or 
arise from Hartree-Fock corrections to the Hamiltonian, 
which have to be computed using nk,f,cr(r, t). 

The Floquet-Boltzmann equati on and the formu¬ 
las for the collision integral (39),p0| m are the main 
results of the first part of the paper. 



FIG. 3. (Color online) Scheme of the hexagonal lattice nsed to 
realize the Haldane model [T]. For every tnnnelling amplitude 
tj {tf) there is an associated lattice vector Vj (uj). Note that 
tf = tf and that the phase of complex hopping strengths is 
defined along the direction of the respective vector in this 
figure. 


V. HALDANE MODEL 
A. Model 

In the following we want to apply the Floquet- 
Boltzmann equation to a concrete example. In a recent 
experiment with ultracold atoms in an optical lattice, the 
Haldane model was realized by means of periodic shaking 
of the lattice [T]. The Haldane model is the prototypical 
example of a topological insulator: Haldane showed that 
an integer quantum Hall state can be realized without any 
external magnetic field on average, but just by arranging 
complex hopping parameters on a hexagonal lattice |36j . 

The experiment can be described (see supplementary 
information of Ref. [1]) by a (distorted) honeycomb lat¬ 
tices, see Fig. with two sites per unit cell, which form 
two chequerboard sublattices A and B. The static Hamil¬ 
tonian can be described by real nearest-neighbour and 
next-nearest-neighbour hopping amplitudes 

H = - ^Uvo.a^u+vc-r) (43) 

U^A,(T 

+ +h-c-) 

+ [Jj' ^Ii+u,, .gQua- + ^iL+vo-(-Uj/,(T^u+vo,o- + h.C.) 

i'.o’ 

with (T being a spin index, and vectors connecting 
points on the same sub-lattice and vectors Vj that con¬ 
nect points on different sub-lattices. 

A periodic shaking of the lattice leads to an accelera¬ 
tion of all atoms. In the frame of reference comoving with 
the lattice, an acceleration can be viewed by a force F(t), 
which is constant in space. Due to the periodic shaking, 
the force is periodic in time, F(t) = F(t -|- Tq). Within 
one period Tg, F(t) rotates in the plane of the lattice on 
an ellipse. The force can either be implemented in the 
Hamiltonian by a potential or, more conveniently, by a 
vector potential A(t) = — (sin(Ht)ei-l-sin(Ht—t/?)e 2 ) 
with dtA{t) = F(t). Here, A is the wave-length of the 


laser used to create the optical lattice, Kq = 0.7778 
parametrizes the strength of shaking in the two perpen¬ 
dicular directions ei and 62 , and H is the oscillation fre¬ 
quency used in the experiment. The parametrization is 
chosen such that A will cancel in the final result. 

The vector potential modifies the hopping amplitudes 
by complex phases. Using Eq. (44), this leads to a de¬ 
scription of the system with complex hopping ampli¬ 
tudes depending periodically on time, Jjlt + To) = Jjit), 
jf^^it + To) = jf^^{t) where 


T _ iZj sin{Qt+^j) j 

Uj — e uj 

jA _ ^izf jA 



(44) 


with 


(A) ^ fKo (A) 

fo A 

= Vj • ei -I- Vj ■ 62 

= Uj • ei - 1 - Uj • 02 


(45) 


and pj > 0. We have implemented the equations for 
two sets of parameters. In appendix 0 we show the 
results for the microscopic Hamiltonian studied in the 
main text of Ref. [1], which is anisotropic and includes 
both nearest-neighbor and next-nearest neighbor inter¬ 
actions. All qualitative features are, however, unchanged 
compared to a simpler set of parameters investigated in 
the following. To mimic the situation discussed in the 
supplementary material of Ref. [T] (see also Ref. [37]), 
where a spinfull, interacting system has been studied, 
we set Jq = Ji = J 2 = J = —2-k ■ 172 Hz, and use 
that approximately Jq = = Aab = 0 and set 

(j) = 7r/2, thus describing a situation where the ground 
state is in the (quantum-Hall) topological phase (note 
that in the experiment also a hopping coupling adja¬ 
cent honeycomb layers was present, which is, however, 
ignored here). Furthermore, we use [T] Vg = A(0.438,0), 
vi = A(-0.062,0.5) and V 2 = A(-0.062,-0.5). We only 
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consider the case of a translationally invariant system at 
half filling. 

To analyze theoretically the effective Hamiltonian aris¬ 
ing from the shaking of the lattice, the authors of [ 1 ] used 
the so-called Magnus expansion, equivalent to a determi¬ 
nation of the Floquet eigentstates to first order pertur¬ 
bation theory in 1/H. Energy non-conserving processes 
can, however, not be treated with a simple Magnus ex¬ 
pansion. 

We will now consider the effects of interactions by 
adding to the Hamiltonian the term 


-ffint = UJ2i (46) 


with Uia = = &l.+vo.a V+vo^ ) ^he A 

(B) sublattice, respectively. While in Ref. [T] mainly the 
non-interacting spinless case was considered, the authors 
also briefly studied the spinfull limit where the interplay 
of interactions and periodic modulations is expected to 
heat up the system. 

The Hartree correction turn our to give only tiny con¬ 
tributions. It leads to a periodically oscillating term 
BAB{t) ~ cU cos(Ht) where c < 0.2 depends on the oc¬ 
cupation function of all states. The Hartree contribution 
therefore remains small compared to all other terms for 
values oiU J where our perturbative formulas can be 
applied. The single-particle gap, for example, changes 
only by 0.2% for 17 = J and c = 0.2. We have therefore 
neglected the Hartree correction which considerably sim¬ 
plifies the numerics as all scattering matrix elements have 
to be computed only once. Furthermore, this approxima¬ 
tion implies that the dependence of U can be absorbed 
into a redefinition of the time, see below. 

To determine the Floquet eigenstates and energies for 
the Floquet-Boltzmann equation it is therefore sufficient 
to diagonalize the non-interacting Floquet Hamiltonian 
given by Eq. (18), where the Fourier components (for 
each spin species) can be written as 2 x 2 matrices (for¬ 
mulated in momentum space) 


/,U _ lU H 

'''n ~ ' '^n,x'^oc 


+ ^n,y^y 


(47) 


where ai are Pauli matrices and 

3 ' 

h^,x = IY (48) 

3 

3 

Here, Jn{') describes the Bessel function of the first kind 
to nth order. Note that the formula for ^ corrects a 
typo in the supplement of [1]. Numerical diagonalization 
of the Floquet Hamiltonian yields then the Floquet states 
4>^ ^ and energies ek,^. 




/Cy(\/27r/a) 


FIG. 4. (Color online). Two-particle collision processes for 
(a) the energy conserving case (n = 0) and (b) the energy vio¬ 
lating case (n = 1). The energy bands of the system described 
in the main text (fl = 6.3| J| = 27r • 1080 Hz) are shown along 
the diagonal of the qnadratic Brillouin zone (fe^, = 0 within 
our conventions). Scattering is indicated by arrows from ini¬ 
tial states (black) into hnal states (blue/red). 


B. Quasi-equilibrium and heating rate 


We consider a translationally invariant situation with¬ 
out external forces. Hence, the Floquet-Boltzmann equa¬ 
tion in Eq. (42) reduces to the simplified form 


— .^coll (49) 

with nk, 5 ,^(t) = rik,{,t(^) = ''T'k,^(i)- The Floquet states 
and quasi energies appearing in the collision integral are 
time-independent. 

To study the heating of the system, we consider the 
change of the (Floquet-) energy par lattice site defined 
by 



dk 

{2TTlaY 


E! ?^k,J,o-(f) 
C.<r 


(50) 


where the factor 1/2 arises as there are two lattice sites 
per unit cell. The heating rate, 7 = dE/dt, is therefore 
given by 


7(t) = 


dE{t) 

dt 


1 

2 


dk 

(27r/a)^ 


hk.c.a(t) (51) 

5,<r 


Using the right-hand side of the Floquet-Boltzmann 
equation (42) and the symmetry properties of the scat¬ 
tering rates, the heating rate can be written as 


_ f ^ <^^2 dqs 

X nQ <^(ek,5 + ~ ~ ^qa.A ~ nfl) 

kqiq2q3 

X ^ (27r/a)^ (5(k -b qi - q2 - qs - crG) 

X 'nq2,riit) nq,^^x(t) (1 — nk,{(t)) (1 — ?T'qi,/2(0) 

( 52 ) 
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FIG. 5. (Color online) Heating rate per lattice site (in units of 
plotted against dimensionless temperature T/ J, where 
J is the hopping amplitude of the isotropic Hubbard model at 
hand. Different curves describe different driving frequencies 
H, ranging from 5| J| = 2?! • 864Hz to 8.8|J| = 27r • 1512Hz, 
see legend. The inset shows a double logarithmic plot. 



FIG. 6. (Color online) Specific heat per lattice site, c(T), 
plotted as function of dimensionless temperature T/J for a 
driving frequency Q. = 6.3| J| = 27r • 1080Hz. The inset shows 
the corresponding entropy per lattice site {kB ~ 1) as function 
of T/J. c{T) depends only very weakly on Q, in the considered 
parameter regime. 


The energy changes in quanta of nfl, determined by the 
Floquet matrix elements IT" and the occupation func¬ 
tions. 

The occupation functions 'Uq 2 , 7 )(i) can be determined 
from the solution of the Floquet-Boltzmann equation. 
This is in general a formidable task due to the high¬ 
dimensional integrals occurring in the collision inte¬ 
gral (39). The problem can, however, be simplified dra¬ 


matically in situations where the scattering rate, 1 /rcon, 
of energy-conserving scattering processes (n = 0 ) domi¬ 
nates over the scattering rate, I/tvIo, for processes which 
violate energy conservation 


-» - 

^con ^vio 


(53) 


In Fig. I^we depict two typical processes associated with 
both time scales. For the parameters studied experimen¬ 
tally in Ref. [T] these rates differ by much more than an 
order of magnitude as discussed below. 

Energy-conservating processes lead to equilibration. 
This implies that under the condition of Eq. (531 after 


a few energy-conserving scattering events the occupation 
functions are well approximated by thermal ones 


nk, 5 (t) « nk.c(T’W) = exp 


m 


+ 1 


-1 


(54) 

where /i(T) is generally determined from the condition 
that the total number of particles remains constant. Note 
also that we set ks = h = 1. 

To leading order in Tvio/Tcon we can therefore replace 
all occupation functions both in Eq. (51) and Eq. (52) by 


Fermi functions, n(t) n°(T(t)). We therefore obtain 


dT{t) dE 


lim) 


(55) 


which can directly be solved by 


t = 


It 7 ( 7 "') 


(56) 


where c{T) = dE/dT is the specific heat. Hence, under 
condition ( |^ we do not have to solve the complicated 
coupled integro-differential equations (H^ for nk,^ (f) • In¬ 
stead, it is sufficient to calculate 7 and E as function 
of temperature and to solve a simple one-parameter dif¬ 


ferential equation (55). Note that the high-dimensional 


integral (52) can still be numerically demanding. 


To determine 7 ( 7 ’), one first has to calculate the matrix 
elements. For the experimental parameters the Floquet 
eigenstates are well localized. For driving frequencies of 
H > 3| J| « 27r-500Hz we have checked that it is sufficient 
to keep track of a few Floquet modes only, e.g., Np = 3. 
Note that IT" with \n\ > 2 cannot contribute to Eq. 


(52) as Ickfl < fI/2. To calculate the integrals in (52) 
taking into account energy conservation, we discretize not 
only momentum by a 20 x 20 mesh but we also discretize 
the energies by rounding them to multiples of AE = 
2tt ■ 12Hz « 0.07|J|. We have checked that for these 
parameters discretization errors are negligible. 

Fig. [^shows the heating rate j{T) for different driving 
frequencies fl. One observes that all heating rates show 
the qualitative same behaviour: 7 starts out at some 
maximal value at T = 0 and then approaches 7 = 0 for 
T —>■ 00 . To analyze the limit T —>■ 00 , it is useful to 


realize that in the situation where (54) holds, a detailed 
balance condition relates the rate T) which is defined 
by the terms in 7 proportional to IT“'"^, to 7 “, the terms 
proportional to IT“^ 


1 _ = gO/T 

7 


(57) 


This follows from the identity for Fermi functions 
n(ln°(l ~ ~ ~ ’^«)(1 ~ n°)e^^ for 


dt dT 
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FIG. 7. (Color online) Dimensionless temperature, TjJ, of 
the interacting system as a function of time for different driv¬ 
ing frequencies Q = 2tt x 864,..., 1512 Hz), see legend of 
Fig. H The inset clarifies that the curves undergo two distinct 
regimes before rising with exponential speed. 



FIG. 8. (Color online) Entropy per lattice site (per particle) 
as a function of time for different driving frequencies 12 (H = 
27r X 864,..., 1512 Hz), see legend of Fig. The solid curve 
(12 = 6.3|J| = 27r • 1080 Hz) is the lowest frequency used in 
the experiment [I]. 


-I- — ca = f2. As 7 = 7 ''' — 7 , we obtain 

liT) oc ^ (58) 

for T » 11 . 

For increasing driving frequency, the heating rate drops 
rapidly, see Fig.[^ This has two reasons: first, the matrix 
elements drop for increasing 12 with oc for 
12 —>■ 00 . Second, the phase space for two-particle scatter¬ 
ing with energy transfer ±12 vanishes due to the restric¬ 
tions on energy conservation for 12 > 211(12). Here, 11(12) 
is the total bandwidth, 11 ( 12 ) = maxk,^ ek{ ~ mink,{ Ck^- 
By expanding around the band minima and maxima, we 
obtain 

7 oc (Hmax - (59) 

for 12 y Umax, ^ flmax with 2Zl(12rnax) — f^max- 
As for 12 > 12max, two particle scattering cannot con¬ 
tribute to heating, one has to consider higher-order scat¬ 
tering events. Within our model, we obtain Umax = 
27r • 1778.5Hz = 10.34 |J|. 

The specific heat per lattice site is determined from 



Fig.§ shows the specific heat as a function of tempera¬ 
ture. For T —>■ 0 the specific heat is exponentially sup¬ 
pressed due to the band gap Aq with Aq « 0.3 J for 
12 = 27r-1080 Hz (for all 12 considered by us the system re¬ 
mains in the gapped topological state). For Aq <T < J, 
the system is approximately described by a Dirac equa¬ 
tion and c(T) grows with T^. For T ^ J, m contrast, 
one obtains c{T) oc 1/T^. 

Using either Eq. ( [^ or Eq. ( [^ , we can directly com¬ 
pute the evolution of temperature as function of time (as¬ 
suming Ti = 0 as the initial temperature). This is shown 
in Fig. for different values of the driving frequency 12. 


Due to the exponential suppression of c(T) for T —>• 0, 
initially T{t) rises very rapidly proportional to 1/ log[l/t] 
for T{t) small compared to the band-gap Ag, followed 
by T{t) oc for Aq < T < J (see inset of Fig. [^. 
For larger temperatures, there is an intermediate regime 
with an approximately linear rise of T. Finally, for T ^ 
12, we obtain dT/dt cx T from 7 cx l/T, c(T) cx 1/T^ 
and Eq. ( |55| ), and therefore an exponential rise of the 
temperature 


T{t) oc for T > 12 (60) 


as confirmed numerically. 


Due to the strong dependence of 7 on the driving 
frequency 12, also T{t) depends strongly on 12. For 
12 —12niax, for example, the prefactor c in Eq. (60) 
descreases rapidly, c oc ( 12 max — U)^. 

In cold atom systems, it is not easy to determine 
the temperature of the system. An observable, which 
is sometimes better accessible is the entropy per par¬ 
ticle (also determined in the extended data section of 
Ref. my The reason is that the entropy remains ap¬ 
proximately constant when all optical lattices and in¬ 
teractions are slowly switched off. We therefore show 
in Fig. also the entropy as a function of time. Due 
to the exponential rise of T with time, Eq. (60), the 
entropy approaches its T —>■ 00 limit with exponential 
speed, S « ln[4] — cojT'^ « ln[4] — for t —>• 00 . 

Finally, we have to check whether in the experimen¬ 
tal system the condition on the scattering rates (53) is 


fulfilled. Taking into account the high-temperatures of 
the experimental system, we will check the condition by 
comparing the ratio of the scattering rates averaged over 
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fi/(27r) 


FIG. 9. The ratio of scattering times, Eq. (621, for energy 
non-conserving (W^^) and energy conserving processes (VF°) 
plotted as a function of the driving frequency fl. For large fre¬ 
quencies, fl > 6| J| « 27r-1000Hz energy conserving scattering 
dominates by several orders of magnitnde. Inset: donble- 
logarithmic plot. 


all bands and momenta in the limit T —^ oo. We define 


lAi" _ f '^‘^2 dq^ 

a (27r/a)2 i2n/ar (2V«)^ (2^^)^ 

^ ~ ^<i2,V ~ ^q3.A ~ nfl) 

kqiq2q3 

X (5(k-b qi - q2 - qs - aG) (61) 


and estimate 


^con 

'^vio 


W+i -bW-i 

wo 


(62) 


In Fig. we show numerical values for the ratio (62) 
as a function of fl. For the (smallest) experimental val¬ 
ues, n « 27r • 1000 Hz, the energy conserving processes 
already dominate by almost two orders of magnitude, 
which justifies the approximation of Eq. ( |54[ ) with high 
precision. As discussed above, all two-particle processes 
violating energy conservation die out with (Umax ~ f2)^ 
for n —Umax; see inset of (62). 


VI. DISCUSSION AND CONCLUSIONS 

In this paper, we have derived the Floquet-Boltzmann 
equation for periodically driven Fermi systems starting 
from the Keldysh-dynamics of Green functions. 

It is instructive to compare the effect of the breaking of 
translational symmetry in time by oscillations periodic in 
time with the effect of the breaking of translational sym¬ 
metry in space by a potential periodic in space. Many 
effects are similar: in the first case energy is conserved 
only modulo hQ, while in the latter case momentum is 
conserved only modulo reciprocal lattice vectors, TiG. 
This leads to a heating of the system to infinite tem¬ 
perature and to the decay of any macroscopic momenta. 


respectively. In both cases the system relaxes to a state 
with maximal entropy allowed by the remaining conser¬ 
vation laws. An important difference is, however, that 
in the case of a periodic potential, one obtains a nomi¬ 
nally infinite number of electronic bands, each of which 
has to be described by an extra quantum number and a 
corresponding semiclassical occupation function. In the 
Floquet case, in contrast, we used the same quantum 
numbers as in the energy conserving case and did not in¬ 
troduce any new occupation function. This can be traced 
back to the very different role taken by space and time 
in single-particle quantum mechanics: while the first one 
is promoted to an operator, this is not the case for time. 

The perturbative Floquet-Boltzmann equation derived 
by us can also be extended to the limit of strong interac¬ 
tions in situations where the number of excitations (e.g., 
doublon and holon excitations in a bosonic Mott insu¬ 
lator) remains small. In this case, the transition rates 
on the right-side of the Boltzmann equation have to be 
computed from the solution of the 2-particle scattering 
problem in the presence of periodic perturbation. 

The absence of energy conservation, ultimately heats 
up closed quantum systems to infinite temperature. Our 
calculations have shown that this effect can be quite 
strong for experimentally relevant parameters [T]. Con¬ 
sider, for example, in Fig. the case = 6.3| J| = 
27r • 1080Hz for moderate interactions, U = \J\. In 
this case, the entropy per site rises from 0.5A:Bln[4] to 
0.75fcBln[4] within At « I8/J « 20ms which is short 
compared to the loading times used in the experiment of 
Ref. [T] . This clearly shows the importance of the heating 
effects. A direct quantitative comparison to the heating 
rates observed in the experiment [1] is, however, not 
possible. In Ref. [T] rather large values of C/ = 10 J (and 
U = 20 J), close to (or in) the Mott insulating phase, were 
investigated, which are far beyond the applicability of our 
approach. Also, a three-dimensional coupling was finite 
in this experiment. For a quantitative description of the 
heating rate in the experimental setup it would also be 
necessary to treat the trapping potential which leads to 
inhomogeneous heating and heat transport through the 
trap. The entropy increase of only a few ks/s reported 
in Ref. [T] appear to be rather small compared to the 
values estimated by us for small initial entropies. We be¬ 
lieve that this can only be explained by the fact that the 
initial entropy per site in the lattice was rather high. In¬ 
deed, in Ref. m , which describes a similar setup without 
shaking, entropies per particle of 1.5 fcs « A:Bln[4] and 
2.5 ks « I.81n[4] have been reported before and after 
loading the system into the trap, respectively. 

As we have shown, the heating rate can efficiently by 
controlled by moderate changes of the driving frequency. 
For n larger than twice the total bandwidth heating 
by two-particle collisions is completely absent and only 
higher-order processes (not included in our analysis) can 
take place which are expected to be strongly suppressed 
by the Pauli principle. Increasing the driving frequency 
implies, however, that all effects of the periodic mod- 
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ulation are also suppressed. For the model considered 
by us, the gap of the topological insulator scales with 
l/n for large Q. Increasing the driving frequency from 
6.3 I J| « 27r • 1080 Hz to ~ 10.3 |J| « 27r • 1800 Hz, 
where 2-particle processes are completely suppressed, re¬ 
duces the gap from 0.3 J to 0.15 J. 

For the design of interacting Floquet systems it will 
be important to control the heating processes. Here we 
hope that the Floquet-Boltzmann equation and variants 
thereof can be a useful tool. 
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Appendix A: Floquet-Moyal expansion 


In the following we sketch the explicit derivation of 
the Floquet-Moyal expansion. Starting point is the time- 
convolution of two two-point functions, given by 


= J dti 


(AI) 


/ 27r 

n,m 

(A2) 

In order to find the explicit form of t) one inserts 

the inverse Floquet-Wigner expressions (see. Eq. (28)) 
for B and C into Eq. (lAll). The formula then reads 


n ,m 
n" .m!‘ 


A(t t') = N ' f ^ —i(uj'+n'Q)t i{uj"+m"Q)t' 

^ ^ 27r 27r 

' ,m" r , , // // 

X / dti 

(A3) 

The first technical challenge enters now via the convolu¬ 
tion in the time-argument ti . To make progress here, one 
expands the functions B and C around ti = t' and ti = t, 
respectively. This leads to the following expression 

4^. _ f doj' duj" V- 

^ ^ ^ md.¥¥ [kj [k' 

k=0k'=0 \ \ 

II n 7 / 

n ,m 

X J + {'rn'— n'')0,)ti ^l+l'—k — k' 

(A4) 


where the identity (H — tY = Ylk=o 
been used, and B^’^') describes the 1th (0th) derivative 
of the function B with respect to time (frequency). We 
can now perform the ti integral by using J = 

(—I)“a!/(ja;i)“ / = 27r(—I)“Q;!/(itLii)“(5(a;i). 

One can further use the identity for derivatives of the 
(5-function, i5("^(a;) = (—I)"n!/a:”5(a;). At the same 
time, we want to move under the frequency integration, 
/ di^, the differentiation away from the (5-function. 
The formula can then be written as 


duj' 


AM=l^^d." E EE 


n',m',l fc=0 fe'=0 
n" ^m" X 


1 I n 

nv\¥W\k 


X {-lY+'^d-k'^k^:k'_ J, 

(O./'-fc') 


X Y'-^' (e-d‘-'+n'u)tQm (w', t)) 
X c'A,;YY„ (w", t) 


(0,l-k) 


(A5) 


Due to the fact that — ^ < w < ^, the argument of 
the ^-function above can only be zero if uj' = uj" and 
m' = n". Hence, the energy integration / dw" YY'L 
be straightforwardly performed. In the following, one 
wants to use the identity (/g)^”^ = YYYk=o 
and simplify summations. Eventually, after some steps a 
relabelling of indices yields the form 




duj 






1 , 1 ' 


(A6) 


According to Eq. (A2) we can read off the explicit form of 
Anm{^,t). After a few more steps of simplification, one 
finally finds an expression for the Floquet-Moyal product 

Anm{ui,t) = 

m' 

(A7) 

where is an operator acting only on object B or C, 

respectively. The Floquet-Moyal expansion is obtained 


by expanding the exponential fucntion in A7 One can 
see that the form of this expression is very similar to the 
ordinary Moyal product. Here, however, the formula re¬ 
quires an additional matrix multiplication in the Floquet 
indices taking care of the fast oscillations. Note again 
that this formula only holds in the limit where the in¬ 
verse Floquet-Wigner transformation (28) is valid. 


Appendix B: Semiclassical dynamics and Berry 
corrections 

In order to describe the semiclassical dynamics of the 
Floquet occupations i.e., the left-hand side of 
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the Floquet-Boltzmann equation, including leading or¬ 
der Berry phase corrections, we follow the derivation by 
Wickles and Belzig [33]. Below we will sketch the idea 
and show why their arguments can be analogously used 
for the Floquet picture. To acquire kinetic equations for 
the individual Floquet states v, one will additionally have 
to decouple the QKE by going into the space that diag¬ 
onalizes the Floquet Hamiltonian. This unitary trans¬ 
formation is, of course, given by the matrix of instanta¬ 
neous Floquet eigenstates. This procedure, however, has 
a peculiarity: one diagonalizes the objects only after the 
Floquet-Wigner transformation. In this form objects are 
allowed to be functions of two canonical variables, which 
is (semi-)classically allowed, but quantum mechanically 
forbidden. So in order to preserve more of the quantum 
mechanical nature of our problem we desire a unitary 
transformation that diagonalizes the object of interest, 
e.g. already on the level of the initial convolu¬ 

tion 


j 

J,/ft(Hz) 

Vj/A 

J/‘/ft(Hz) 

Uj/A 

0 

-27r • 746 

(0.438,0) 

- 

- 

1 

-27r • 527 

(-0.062,0.5) 

271 • 14 

(0.5,-0.5) 

2 

-27r • 527 

(-0.062,-0.5) 

277- 14 

(-0.5,-0.5) 

3 

-27r • 126 

(-0.562,0) 

277-61 

(0,1.0) 


TABLE I. Summary of all hopping strengths and associated 
lattice vectors as used in Appendix]^ as well as in [Tj. While 
Jj describes the hopping between nearest neighbours with 
■Vj connecting these points (from sublattice A to B), jf de¬ 
scribes the hopping between next-nearest neighbours with Uj 
connecting points on the same sublattice. Note also that 
ei = (1,0) and 02 = (0,1), and that jd = J?. 

can observe that the expression above is nothing but a 
Taylor expansion to first order. Hence, after introducing 
band projected kinetic variables X = x — and 11 = 
7 r -|- hAx ^, one can write 


[/o(G^)-io[/1'= (G^)-i (Bl) 

with UoW = 1 and where ’o’ can be viewed as a Floquet- 
Moyal product. The obvious problem is that U can gen¬ 
erally not be found. However, one can systematically 
calculate corrections order by order. Hence, one writes 
U = C/o(l + Gi + • • •) with Go being the unitary ma¬ 
trix that diagonalizes the instantaneous Floquet Hamil¬ 
tonian, UqH^Uq = H^. Expanding the Eloquet-Moyal 
product yields terms like 

Uod.U^o =9r- iA^ (B2) 

where Berry connections have been defined 

A, = tUoid^U^o) (B3) 


Here, di runs over all possible derivatives (space, momen¬ 
tum, time, energy). The crucial point to realise here is 
that Ai is a matrix in Floquet space. While in the deriva¬ 
tion by Wickles and Belzig this matrix only knew about 
the original bands of the Hamiltonian, here the object is 
fully aware of the underlying Floquet structure. The rea¬ 
son why the Floquet notion is implemented here in such 
a straight forward manner lies in the fact that it was 
possible to reduce the fast oscillations to a simple matrix 
structure. Hence, the procedure of diagonalization by 
means of U naturally introduces the Floquet character. 

Performing the unitary Floquet-Moyal transformation 
(Bl I, keeping only terms up to first order in TiA, and pro¬ 
jecting the expression onto the eigenstates of the prob¬ 
lem, i.e., here the Floquet states, the object (G^)“^ can 
be approximated as 


(G«)-i«(G«)o-i(X,n) (B5) 

So we can indeed express the object in its trivially di¬ 
agonalized form. However, the prize we need to pay is 
the change from canonical to kinetic variables. Thus, 
the derivatives appearing in the Moyal product have to 
be changed accordingly, and the new expression reads 

{n = i) 


o ^ exp[-^(a|ag - agag) -b ^a^H^,,a^] (b6) 

for a product of i? o G where the partial derivative 
a-®/*" acts only on object B or G, respectively. Here, 
/r = 0, • • • ,6 run over all time, position and momen¬ 
tum indices. The Floquet-Moyal product ’o’ also requires 
a multiplication of Floquet matrices B and G. Hence, the 
entries of the object il, the Berry-curvature tensor, are 
in fact matrices in Floquet space. The explicit form (for 
non-crossing Floquet bands) reads 


- 


aM! 


id) 


dA 


id) 


dfi 


dv 


(B7) 


Using Eq. (B3) and the fact that Uq is simply the matrix 


of Eloquet eigenstates one can write down the Berry- 
curvature tensor for the Eloquet state ^ 

(B8) 

i,n 


(G«)-1 =(G«)o' - (B4) 

+ ^{A^,di^HG%^} + 0{{hAf) 

where x = (t, r), tt = (w, —p), A^^ = '^■ViAVi, with 
Vi being a projector, and (G^)g ^ = Go(G^)“^Gg. One 


Note that the projection onto on-shell processes per¬ 
formed on the right hand side of the Boltzmann equation 
is al read y encoded in the procedure above by virtue of 
Eq. (Bl|. At the same time, the instruction to go from 
canonical to kinetic variables is obsolete for the right 
hand side, since there we allow only zero order terms 
to contribute (for which X = a: and 11 = tt). 
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FIG. 10. (Color online) Heating rate per lattice site (in 
units of l/H^) plotted against dimensionless temperature 
r/Jmax, where Jmax is the maximal hopping amplitude of 
the anisotropic Hubbard model at hand. Different curves de¬ 
scribe different driving frequencies H ranging from 4| Jmaxl ~ 
27r • 3000 Hz to 8| Jmaxl ~ Stt • 6000 Hz, see legend. The in¬ 
set shows a double logarithmic plot. The solid curve (H = 
27r • 1080 Hz « 5.4|Jmax|) is the frequency used in the main 
text of [I]. 



FIG. 11. (Color online) Specific heat per lattice site, c{T), 
plotted as function of dimensionless temperature T/ Jmax for 
a driving frequency H = 27r • 4000Hz ~ 5.4|Jmax|. The inset 
shows the corresponding entropy per lattice site (fes = 1) as 
function of T/Jmax. 


Appendix C: Alternative set of parameters for the 
interacting Haldane model 

In Sec. |V]we presented the predictions of our Floquet- 
Boltzmann equation]^ for the recently realized Haldane 
model [1] for an isotropic setup, characterized by a single 
nearest-neighbour hopping amplitude J only. 

Here we repeat the analysis for the (spinfull) model de¬ 
scribed in the main text of [T]. The geometry of the model 
remains the same as described in Sec. |V| see Fig. but 
hoppings are anisotropic and also next-nearest neighbour 
hopping amplitudes are taken into account. A sum¬ 


mary of the experimental parameters of Ref. [T] which we 


also used in our calculations can be found in Tab. HI 



FIG. 12. (Color online) Dimensionless temperature, T/Jmax, 
of the interacting system as a function of time for different 
driving frequencies H (H = 27r x 3000,..., 6000 Hz), see legend 
of Fig. |10[ The inset clarifies that the curves undergo two 
distinct regimes before rising with exponential speed. 



FIG. 13. (Color online) Entropy per lattice site (per particle) 
as a function of time for different driving frequencies H (H = 


27r X 3000,..., 6000 Hz), see legend of Fig. 10 


We then apply the exact same arguments and proce¬ 
dures as in Sec. |V] Note that all approximations and 
assumptions made above also hold here because the driv¬ 
ing frequencies is here now tuned to higher absolute val¬ 
ues n = 2Tr ■ 3000,..., 6000 Hz. Translating these fre- 
quncies into units of the maximal tunnelling amplitude 
Anax = Jo, one sees that the different values of H span the 
same regime as above, i.e., H « 41 Jmaxl, • • ■, 81 Jmaxl- We 
now recalculate Figs. [5@for the alternative set of param¬ 
eters presented here and plot the results in Figs |l0p^ 
Note that we used Jmax throughout to rescale energies. 
In fact, comparing the two sets of parameters clearly re¬ 
veals that the results are qualitatively unaffected. Hence, 
all conclusions drawn from the results in Sec. |V] also ap¬ 
ply here. 
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